clear all;
%% Code explanation
% 170305 DZ 
% 190709 DZ updated to calculate 50%Rn and 1%Rn
% 190712 DZ use 3sigma to estimate Tmc_err
% 200121 DZ corrected the resistance factors, export to spreadsheet data

%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter

L_W1=43/66; % 3Sn-13PbTe
L_W2=95/61; % 5Sn-15PbTe

ratio=0.5;
%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

root=('D:\Documents2\Rawdata\Sn\180801 MPI mK\stanene_3_5_layer_outofplane\m'); %

filenames={'0uW';'20uW';'40uW';'60uW';'80uW';'100uW';'130uW';'160uW';...
    '220uW';'280uW';'340uW';'400uW';'480uW';'580uW';'700uW';'850uW';...
    '1000uW';'1200uW';'1500uW';'2000uW';'2600uW';'3200uW';'4000uW';'4800uW';...
    '4200uW_rot';'5000uW_rot';'5500uW_rot';'6000uW_rot';'7000uW_rot';'8000uW_rot'};



Boffset=0.0127;
NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT,1);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT,2);
Rc=zeros(NT,2);


for kk=1:NT

    file1=[root,filenames{kk},'.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f %f %f','commentStyle','#');
        data1=data1';
     
        fclose(fid);
   
          B=data1{1}-Boffset;
          R1=(data1{2}/L_W1+0.025)/10;
          R2=(data1{3}/L_W2+0.002)/10;
          T=T_RuOx(data1{5}*100);
          Tmc(kk)=mean(T);
          Tmc_err(kk)=3*std(T);
          
          subplot(1,2,1),plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
          hold on
           subplot(1,2,2),plot(B,R2,'Color',colorset(kk,:)','LineWidth',1);
          hold on
       
          if max(B)<2
              BB=1.6;
          else 
              BB=2;
          end
          indexB=find(B>BB);
          
          p1=polyfit(B(indexB),R1(indexB),1);
          p2=polyfit(B(indexB),R2(indexB),1);
          
          index1=find((R1-ratio*p1(1)*B-ratio*p1(2)<=0)); 
        if isempty(index1)==0
          Hc2(kk,1)=B(index1(1));
          Rc(kk,1)=R1(index1(1));
        end
       index2=find((R2-ratio*p2(1)*B-ratio*p2(2)<=0)); 
        if isempty(index2)==0
          Hc2(kk,2)=B(index2(1));
          Rc(kk,2)=R2(index2(1));
        end
         title={'B(T)','R(kOhm)','T(K)','averaged T(K)'};
        AA=[B R2 T Tmc(kk)*ones(length(B),1)];
       
         xlswrite('FigS5A.xlsx',AA,kk);
       xlswrite('FigS5A.xlsx',title,kk);
      
%     
end

p11=polyfit(Hc2(1:24,1),Rc(1:24,1),1);
p22=polyfit(Hc2(:,2),Rc(:,2),1);
subplot(1,2,1),plot(B,p11(1)*B+p11(2),'-k');
              hold on
subplot(1,2,2),plot(B,p22(1)*B+p22(2),'-k');
              hold on

pp1=polyfit(Tmc(19:24),Hc2(19:24,1),1);
pp2=polyfit(Tmc(21:30),Hc2(21:30,2),1);

Tmc1=abs(pp1(2)/pp1(1));
Tmc2=abs(pp2(2)/pp2(1));

xi=sqrt(6.626E-34/4/pi/1.6E-19/(abs(pp2(2))))*1E9;

Tp=(0:0.002:Tmc2)';

figure
plot([Tmc' Tmc1]',[Hc2(:,1)' 0]','dr',[Tmc' Tmc2]',[Hc2(:,2)' 0]','sb',Tp,pp2(1)*Tp+pp2(2),'-k');
hold on





%% End




